Spherical Collapse and the Halo Model in Braneworld Gravity 
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We present a detailed study of the collapse of a spherical perturbation in DGP braneworld gravity 
for the purpose of modeling simulation results for the halo mass function, bias and matter power 
spectrum. The presence of evolving modifications to the gravitational force in form of the scalar 
brane-bending mode lead to qualitative differences to the collapse in ordinary gravity. In particular, 
differences in the energetics of the collapse necessitate a new, generalized method for defining the 
virial radius which does not rely on strict energy conservation. These differences and techniques 
apply to smooth dark energy models with w ^ — 1 as well. We also discuss the impact of the 
exterior of the perturbation on collapse quantities due to the lack of a Birkhoff theorem in DGP. 
The resulting predictions for the mass function, halo bias and power spectrum are in good overall 
agreement with DGP N-body simulations on both the self-accelerating and normal branch. In 
particular, the impact of the Vainshtein mechanism as measured in the full simulations is matched 
well. The model and techniques introduced here can serve as practical tools for placing consistent 
constraints on braneworld models using observations of large scale structure. 



I. INTRODUCTION 

Modified gravity models have attracted a great deal 
of interest recently as an alternative explanation of the 
observed accelerated expansion of the Universe [iJ-IH . In 
order for this scenario to work, gravity must be signifi- 
cantly modified from General Relativity (GR) on cosmo- 
logical scales, but has to reduce to GR locally in order to 
satisfy stringent Solar System constraints at a few AU. 
Thus, a working modified gravity model has to include a 
non-linear mechanism to restore GR in high-density en- 
vironments, which can have a noticeable impact on the 
formation of large-scale structure on intermediate scales 
of a few to tens of Mpc . 

One popular modified gravity scenario is the DGP 
braneworld model (lpj . Here, a four-dimensional 
Friedmann-Robertson- Walker universe is imbedded as 
a brane in five-dimensional Minkowski space. In this 
model, gravity is 5-dimensional on the largest scales, and 
becomes four-dimensional at scales below the crossover 
scale r c , a fundamental parameter of the model. The 
modification of the Friedmann equation depends on the 
choice of embedding for the brane [ll| which yields two 
branches of the model: the self- accelerating branch, on 
which accelerated expansion occurs at late times with- 
out any cosmological constant or dark energy, and the 
normal branch where no such acceleration occurs, and a 
brane tension or other form of stress-energy with negative 
pressure has to be added on the brane. 

On scales much smaller than the horizon and crossover 
scales, DGP gravity can be described by an effective 
scalar-tensor theory [l2|, EH , where the additional scalar 
degree of freedom, the brane-bending mode ip is asso- 
ciated with displacements of the brane from its back- 
ground position. The brane-bending mode yields an ad- 
ditional gravitational force which influences dynamics of 



non-relativistic particles. On the normal branch of DGP 
this force is attractive, while it is repulsive in the self- 
accelerating branch. Non-linear interactions of p, via 
the so-called Vainshtein mechanism, ensure that it has 
tiny effects within the Solar System. For values of r c of 
order the current horizon, these interactions become im- 
portant as soon as the density contrast becomes of order 
unity. Hence, it is crucial to consistently follow the full 
brane-bending mode interactions in a cosmological simu- 
lation, as has recently been done d, E3 • While we study 
the behavior of <p in specific DGP models, the form of 
the non-linear interactions is expected to be generic to 
braneworld models with large extra dimensions [TBI . [l6j . 

Given the considerable computational expense of these 
simulations, it is worthwhile to develop a model which 
captures the main modified gravity effects, enabling fore- 
casts and constraints properly marginalized over cosmo- 
logical parameters. In the halo model of large structure 
)l7l | , which assumes that all matter is in bound dark mat- 
ter halos, the abundance and clustering of halos are deter- 
mined by the linear power spectrum once characteristic 
quantities of the collapse of a spherical perturbation are 
determined, namely, the linear collapse threshold and the 
virial radius or overdensity. 

Spherical collapse is well studied for General Relativity 
with a cosmological constant (e.g., [3 Ell), an d has also 
been explored for quintessence-type dark energy [20l - [22| . 
A first study of spherical perturbations in the context of 
DGP was undertaken in [6| . We extend previous studies 
by deriving the full ip field profile of an isolated mass in 
closed form, and by carefully defining the interior and 
exterior forces during collapse along with the energetics 
implied by the profile. In particular, the potential energy 
required for the virial theorem and for the total Newto- 
nian energy differ, with the latter not being strictly con- 
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served during collapse. 1 This lack of a conservation law 
also applies to dark energy models with an equation of 
state w — 1. We show how this problem can be cir- 
cumvented by properly defining the condition for virial 
equilibrium based on forces. 

We discuss the general properties and parameteriza- 
tion of DGP models in § HH and the spherical collapse 
calculation in § IIIII The halo model calculations are out- 
lined in § IIV1 and the results and comparisons with sim- 
ulations are given in § |V] We conclude in § I VII The 
Appendix contains derivations of the tp profile and other 
quantities needed in the collapse calculation, as well as a 
discussion of the potential energy and virial theorem in 
DGP. 



II. DGP MODELS 

In this section, we discuss the general properties of 
the DGP models considered in this paper. The first 
model (sDGP) is in the self-accelerating branch of DGP 
with neither a cosmological constant nor spatial curva- 
ture. During matter domination and beyond, the modi- 
fied Fricdmann equation in sDGP reads: 



H 



sDGP 



(a) = H 



fi. 



where 
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(1) 



(2) 



and po is the average matter density today. This expan- 
sion history is clearly different from ACDM, and corre- 
sponds to an effective dark energy with w e g — > —1/2 in 
the matter-dominated era at high redshifts. For com- 
parison, we will also consider an effective smooth dark 
energy model (QCDM) with the same expansion history 
as sDGP. Note that this expansion history when com- 
bined with the growth of structure near the horizon is 
in substantial conflict with data [23| . Moreover, the self- 
accelerating branch is plagued by ghost issues [H, HH, [25| 
when perturbed around the de Sitter limit. Despite these 
problems, sDGP remains an interesting toy model for ac- 
celeration from modified gravity. 

The second scenario {nDGP) is in the normal branch 
of DGP. In order to achieve acceleration, it is necessary to 
add a stress-energy component with negative pressure on 
the brane. We adopt the model introduced in [26| . where 
a general dark energy component is added on the brane 
but the geometry remains spatially flat. The equation of 
state of this dark energy is adjusted so that the expansion 



1 The "total energy" here is defined for a Newtonian cosmology 
and its non-conservation does not imply violation of covariant 
energy-momentum conservation. 



history is precisely ACDM: 

H nDGP (a) = Ho^/flma- 3 + f2 A , (3) 

again during matter domination and beyond. While fl m 
quantifies the true matter content in this model, is 
to be seen as an effective cosmological constant relevant 
for the expansion history only. This construction allows 
our nDGP models to evade the otherwise stringent ex- 
pansion history constraints on r c with a true cosmolog- 
ical constant or brane tension (2~7j . Likewise it provides 
a class of models where the observable impact of force 
modification is cleanly separated from the background 
geometry. Wc consider the two models of [26}, nDGP-1 
with r c = 500 Mpc and nDGP~2 with r c = 3000 Mpc. 

For notational simplicity, it is convenient to also phrase 
the sDGP Fricdmann equation in terms of an effective 
dark energy component so that in both cases 
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(p + Pes), 



(4) 



where p is the background matter density and p e s is im- 
plicitly defined by Eq. flU and Eq. ©. In Tab. [fl wc 
summarize the parameter choices for the simulated mod- 
els 0, Hi| ; in case of sDGP, they are from the best-fitting 
fiat self-accelerating model of [H to WMAP 5yr, Su- 
pernova and Ho data, while for the nDGP models, the 
expansion history and primordial normalization match 
those of the best-fitting ACDM model of 0]. 

On scales much smaller than both the horizon H~ x , 
and the cross-over scale r c , DGP reduces to an effective 
scalar-tensor theory with the brane-bending mode <p rep- 
resenting the scalar. Time variation in ip induced by the 
nonrclativistic motion of the matter involve the dynam- 
ical time and can be neglected with respect to spatial 
derivatives in this regime. The tp field then couples to 
matter by contributing to the metric potentials ^, <I> de- 
fined by the line element 



(1 + 2^)dt 2 



a 2 (l 



2<S>)dx 2 



as 



(5) 

(0) 
(7) 



where ^jv is the Newtonian potential determined via the 
usual Poisson equation 



V 2 f N =4irG5p, . 



(8) 



Here and throughout, spatial derivatives are physical, not 
comoving. 

While the motion of massive, non-relativistic particles 
such as cold dark matter is governed by the dynamical 
potential ^P, the propagation of light is determined by 
the lensing potential — $)/2. This combination is not 
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affected by ip due to the conformal invariance of electro- 
magnetism. Hence, in DGP lensing mass is equal to the 
"actual" mass, while the dynamical mass differs unless 
|<p/tfjv| < I- 

In the sub-horizon, quasi-static regime, the equation 
for the brane-bending mode can be written as (e.g., [28j |) 



vV 



3/3 



[( VV) 2 - (V,V^)(V ! VV)] 



8ttG 
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where the function /3(a) is given by 



(a) = 1 ± 2#(a) r c 1 



g(°) 
3# 2 (a) 



<y Pl (9) 



(10) 



Here, the positive sign holds for the normal branch 
of DGP, while the negative sign holds for the self- 
accelerating branch. Note that the sign of /3 determines 
whether the force mediated by the brane-bending mode 
is attractive {fi > 0, nDGP) or repulsive (/3 < 0, sDGP). 



so we can ignore the impact of any compensating undcr- 
density on the dynamics of collapse. Note that this is 
unique to a top-hat density and does not hold for other 
density profiles, in which case the collapse will not be 
self-similar anymore. 

Given p and R, there are two important mass parame- 
ters: the total mass M = 4npR 3 /3 and the mass pertur- 
bation SM = AirSpR 3 /3. The first is conserved during 
collapse, while the second is the source of the <p field 
and gravitational potential, and is therefore useful in ex- 
pressions involving the field profile. With the definition 
S = 5p/~p, the two masses are related by 



SM 



1 + 6 



(12) 



such that they coincide at high overdensity. 



A. Collapse Dynamics 



III. TOP-HAT COLLAPSE IN DGP 

In this section, we review the dynamics of the collapse 
of a top-hat density perturbations in the DGP case. We 
follow 0, [§[ but pay special attention to the ip profile as 
well as subleties in the potential energy and virial condi- 
tion which arc specific to the Vainshtein mechanism (see 
Appendix for details). 

We assume an initial top- hat density profile of the form 



p(r) 



Sp, r < R, 
0, r>R, 



(11) 



We show in the Appendix that the top-hat profile remains 
top-hat during the collapse despite sweeping out an un- 
derdensity outside of R. We further show that forces in- 
side of R depend on the enclosed mass perturbation and 



TABLE I: Parameters of the simulated DGP cosmologies. 





QCDM 


sDGP 


ACDM 


nDGP-1 


nDGP-2 




0.258 


0.258 


0.259 


0.259 


0.259 


(eff.) 








0.741 


0.741 


0.741 


r c [Mpc] 


oo 


6118 


oo 


500 


3000 







0.138 





17.5 


0.487 


Ho [km/s/Mpc] 


66.0 


66.0 


71.6 


71.6 


71.6 


too n b h 2 


2.37 




2.26 




ft c h 2 


0.089 




0.110 




r 


0.0954 




0.0825 




n s 


0.998 




0.959 




,4 a (0.05 Mpc -1 ) 


2.016 10~ 9 




2.107 10" 


9 


cr 8 (ACDM) a 


0.6566 


0.7892 



"Linear power spectrum normalization today of a ACDM model 
with the same primordial normalization. 



Given a metric specified by 'I', conservation of 
energy-momentum is unchanged in DGP and leads to 
the same equation of motion for the density perturba- 
tion as in ordinary gravity. On scales much smaller than 
the horizon 



s-i- 



3 1 + 5 



2H5 = (1 + <5)V 2 *, 



(13) 



where H = a/a denotes the Hubble rate, and dots de- 
note derivatives with respect to time. The modification 
of gravity enters through the dynamical potential ty: in 
general relativity (GR), * = $ N [Eq. ©], while in DGP 
'J receives an additional contribution from the brane- 
bending mode ip following Eq. ([6]). 

The full profile of ip around a top-hat perturbation is 
derived in the Appendix. The key result is that in the in- 
terior of the top-hat, V 2 ip is constant like V 2 "!^. Hence, 
a pure top-hat will stay top-hat, so that the Eq. ((To) can 
be considered as an ordinary differential equation involv- 
ing a spatially constant S. Note that as shown in § IA 71 
this is unique to a top-hat and will be violated as soon as 
more general spherically symmetric profiles are consid- 
ered. Moreover the implied scaling with the local matter 
density is not true for the exterior of the top- hat (cf. [29l|). 
We discuss this issue further in the Appendix. 

While W 2 ip is spatially constant for r < R. it has a 
nontrivial dependence on Sp (see Appendix). This depen- 
dence can be cast in terms of an effective gravitational 
constant: 



V 2 ip = 8^AG DGP (i?/i?,)5p, 

AGdgpW = — — G > 



3/3 



where the Vainshtein radius 
i?* = 



fl&GSMr 2 ^ 1/3 



V 9/3 2 



(14) 
(15) 

(16) 
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If R ^ i?*, AGdgp = G/(3/3) and (p is simply propor- 
tional to the Newtonian potential, ^jv( r ) = G8M/r. We 
call this the linearized limit as it applies to the 8M <C M 
limit. 

In the opposite limit relevant for 5p/~p^> 1, AGdgp <x 
(i?/i?*) 3 / 2 1. This is the Vainshtein limit and here 
ip ~ G5M/R* is to leading order constant throughout 
the perturbation (§ IA 2[) . 

We can use mass conservation 



S) 



(17) 



to rewrite the top-hat equation of motion in Eq. (|13[) as 



- = H 2 



H - -V 2 * 
3 



4ttG 



[p + (1 + 3w e ff )/Oeff] 



47rG D Gp(i?/i?*) 



(18) 



where Gdgp = G + AGdgp- In the second line, the 
effective equation of state is defined to be —3(1 + u> e ff) = 
dlnpoff/dlna. These terms account for the effect of the 
background expansion. 

The terms involving the expansion history can be seen 
as coming from an effective potential obtained by expand- 
ing the Friedmann-Robertson- Walker metric around the 
center of the perturbation [lj| [2l[ : 



cir 



2ttG 



(p+ (1 + 3w e ff)Peff) r 2 



(19) 

up to a constant that is irrelevant for the dynamics (see 
I A 7j) . Note that V 2 5 , e ff is spatially constant and so its 
effect also preserves the top-hat profile. 

Eq. (|18[1 can then be written in compact form as 



| = -i(V 2 vI/ cff + V 2 vI/). 



(20) 



Note that the pieces involving the matter density com- 
bine as 



V 2 (* cff + tfjv) = 4nG[p + (1 + 3w ce )p eS } 



(21) 



so as to reflect the total matter density p inside the top- 
hat as one would expect from Newtonian mechanics (see 
Appendix for further discussion). 



B. Collapse calculation 

We numerically solve the spherical collapse equa- 
tion (fT5|) following Q . Specifically, we start at an initial 
scale factor a, = 10~ 5 , using In a as a time variable, and 
replacing R with y defined as 



V 



R 



a 
a, 



(22) 



where Ri is the initial radius of the perturbation. Hence, 
we start with y = and y' = —Si/3 as given by linear 
theory in the matter dominated epoch in terms of the 
initial density fluctuation Si. Here, we have set AGdgp = 
at dj, since the effects of force modifications in DGP 
are negligible at such an early time. 

With these initial conditions, we can solve Eq. (|18[) as 



H' 

ir 3 



V 



H' 

Gbgp(R/ R*) 



2H 2 {a) 



G 



a 

a, 



S. (23) 



The overdensity relative to the background 5 is given by 



5(y,a) = (l + 6i) (^V + l)" 



- 1. 



(24) 



Fig. [T] shows the result of solving this equation in the 
different models (bottom panel). We adjust Si so that 
collapse, where R = or y = —a/ai, occurs at a = 1. 
Turnaround, where dR/dhia = or y' = —a/ai, occurs 
at a — 0.54 — 0.56 for these models. 

Fig. Q] also shows the evolution of the gravitational 
force strength Gdgp/G. Since collapse is defined by 
S — > oo at a — > ao, the perturbation eventually becomes 
much smaller than its Vainshtein radius in all models, so 
that Gdgp — > G. Thus, the evolution of forces between 
turnaround and collapse is significant. This evolution 
raises the issue of the conservation of total energy of the 
perturbation during collapse. We will return to this ques- 
tion in the Appendix. 

We then extrapolate the initial overdensity Si to ao = 1 
using the linear growth equation, obtained from lineariz- 
ing Eq. CE3]) : 



2H 5 = AnGnnpS, 



where 



Glin(o) 



30(a) 



G 



(25) 



(26) 



is the linearized value of Gdgp for R/R* » 1 in Eq. (|23|) . 
The resulting overdensity is the linearly extrapolated col- 
lapse overdensity, which we call S c . 

To expose the impact of the Vainshtein mechanism, 
we will also consider linearized DGP collapse (note that 
this it not linearized collapse) as the limit of no non- 
linear ip interactions. In this case, we replace Gdgp with 
Gn n in Eq. (f2"3"|) . Note that the Vainshtein mechanism 
is strongest for a spherically symmetric collapse and ab- 
sent for a planar collapse and so these two cases should 
encompass the range of possibilities in the cosmological 
context HIH. 

For a top-hat, the Vainshtein radius in units of R is 
given by = (eS)^ 1 ^ 3 , where e is defined as 



9/3 2 (a) 



(H r c ) 2 VL m a" 



(27) 



5 




a 




a 



FIG. 1: Evolution of gravitational force Gdgp /G as a function 
of a for perturbations collapsing at ao = 1 in the different 
DGP models (top panel). In each case, the thin lines show 
the linearized force modification Eq. (|26[) whereas the thick 
lines show the nonlinear case. We also show the evolution 
of the scaled radius R/Ri of the perturbation in the bottom 
panel. 



such that 1/e represents the density threshold at which 
the top-hat perturbation crosses into its own Vainshtcin 
radius. Fig. [2] shows e and e<5 as a function of a for a 
range between turnaround (a ~ 0.5) and collapse (a = 1) 
for the simulated models of Tab. |U In the nDGP models 
this threshold increases substantially toward the present 
and the Vainshtcin suppression does not saturate until 
quite late in the collapse. In fact in the nDGP-1 model, 
we shall see that the perturbation virializes before the 
Vainshtein mechanism can operate. 



C. Virial Radius and Overdensity 

Spherical top-hat collapse formally predicts a collapse 
to a singularity at R — >• 0. In reality we expect that pro- 
cesses such as violent relaxation will eventually establish 
virial equilibrium. More specifically, the virial theorem 
relates the kinetic energy of the body 

T = Jd 3 x l -pv 2 = ^MR 2 , (28) 

where the last equality holds for a top-hat, to the trace 
of the potential energy tensor 

2T + IV = 0, (29) 



FIG. 2: Top panel: The non-linearity parameter e. The Vain- 
shtein mechanism operates once e S > 1. Bottom panel: e S for 
a top-hat perturbation that collapes at ao = 1 in each model. 



with the definition 

W ee - J Ap m (x)x-Vf. (30) 

Note that W depends explicitly on forces only and some 
care must be taken in relating it to the potential or bind- 
ing energy of the perturbation for contributions from the 
brane bending mode <p (see Appendix IA 3[) . We show 
there that W can be broken up as a sum of three contri- 
butions to the dynamics, 

W = -^-^(l + 3» cff WMfl 2 
5 R 5 

3 AGdgp MSM 

5 R [ ' 

We determine the virial radius i? v ir of the perturba- 
tion as the radius during collapse (after turnaround) at 
which Eq. ([29)) is satisfied. In the literature conserva- 
tion of total energy is often used to set this condition. 
One can easily show that if energy conservation holds 
strictly, our determination of i? v ; r agrees with the usual 
definition. However, in the presence of an evolving p e g 
and AGdgPi energy is no longer strictly conserved over 
a Hubble time. Note that this problem occurs in dark 
energy models as well but is usually ignored under the 
assumption that w c s ~ — 1. However, this procedure is 
not justified when considering modifications due to a fi- 
nite 1 + w e ff. Differences between our definition of i? V ir 
and the one relying on energy conservation are of the 
same order as the differences induced by 1 + iy e ff- We 
discuss these issues further in Appendix I A 5lTA~6l 
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Once i? v i r is calculated, we have the density of the 
perturbation at this radius, p vu = p(a v ; r )[l + (5(i? v ; r )], 
where a v ir is the scale factor at which the perturbation 
reaches i? V j r during collapse. For collapse at do = 1, 
a V j r rj 0.91 — 0.93 for the simulated models. One then 
assumes that the perturbation maintains this "virial den- 
sity" while the background continues to decrease. The 
final virial density with respect to the background, A v j r 
at ao is then given by referring this density to the back- 
ground matter density at oq: 



ir = [1 + S(R vil )} [ 



(32) 



Tab, im shows the resulting spherical collapse parameters, 
S c and A v j r , for the different models and gravitational 
modifications: unmodified (GR collapse), valid for GR 
with smooth dark energy; the expression Eq. (fT5j) (DGP 
collapse)] and Eq. (|26p {linearized DGP collapse). 



IV. HALO MODEL PREDICTIONS 

We now briefly describe how we move from spherical 
collapse predictions of the linear collapse threshold 8 C 
and the virial overdensity A v i r summarized in Tab. |n] 
to predictions of the halo mass function, bias, and non- 
linear power spectrum. For further details, see 

In the Press-Schechter approach, one assumes that all 
regions with 5 > 5 C in the linearly extrapolated initial 
density field collapse to form bound structures (halos). 
The fraction of mass within halos at a given mass is 
then determined by the variance of the linear density 
field smoothed at that scale. Here, we adopt the Shcth- 
Tormen (ST) prescription j3l| for the halo mass function 
predictions, which enables a direct use of our spherical 
collapse results. Also, we previously found a good match 
to the ST mass function and bias in our ACDM simula- 
tions n. 

The ST description for the comoving number density 
of halos per logarithmic interval in the virial mass M v -„ 
is given by 



dn 



dlnM v 



M v 



7M 



dv 



dlnAfv 



(33) 



TABLE II: Spherical collapse parameters for the cosmologies 
defined in Tab. [T] for collapse at ao — 1. 





Collapse type/Model: 


sDGP 


nDGP-1 


nDGP-2 




GR 


1.662 


1.674 


1.674 




DGP 


1.627 


1.687 


1.688 




DGP lin. 


1.676 


1.678 


1.672 


A vir 


GR 


399.9 


372.3 


372.3 




DGP 


467.1 


300.4 


322.8 




DGP lin. 


436.4 


311.7 


339.1 



where the peak threshold v = 5 C J a(M v - a ) and 



vj(v) = A^av 2 [\ + (av 2 )- p ] exp[-ai/ 2 /2] 



(34) 



Here, the virial mass is defined as the mass enclosed at 
the virial radius i? v ir- o~(M) is the variance of the linear 
density field convolved with a top hat of radius R that 
encloses M = 4irR 3 ~p/3 at the background density 



a 2 (R) 



d 3 k 
(2tt) 3 



\W(kR)\ 2 P h (k), 



(35) 



where Ph{k) is the linear power spectrum and W is the 
Fourier transform of the top hat window. The nor- 
malization constant A in Eq. (jMJ) is chosen such that 
j dvf(v) = 1. We adopt the standard parameter values 
of p = 0.3 and a = 0.75 throughout. 

The linear bias corresponding to the ST mass function, 
obtained in the peak-background split, is given by j3l| 



&lin(M v 



= 1 



b(k = 0, M vil 
av 2 - 1 



2p 



S c [l + (ais 2 )P] 



(36) 



By assuming a specific form of halo density profiles, 
we can rescale mass definitions from the virial mass M v ; r 
to .M200, the mass definition used in the simulation mea- 
surements, as outlined in (32| (again, all overdensities are 
referred to the background matter density). We use this 
approach to compare the scaling relation predictions to 
the simulations in SfV] For the halo profiles, we take an 
NFW form [H, 



PnfwM 



r/r s (l + r/r s ) 2 ' 



(37) 



where r s is the scale radius of the halo and the normaliza- 
tion p s is given by the virial mass Af V j r . We parametrize 
r s via the concentration c v ; r = i? v j r /r s given by [34| 



c vir (M vir ,z = 0) =9 



/M v 



-0.13 



(38) 



where M* is defined via <r(M») = 5 C . Since generally 
Rvir,R200 3> r s , the precise form of the concentration 
relation has a negligible impact on the mass rescaling. In 
the following, when no specific overdensity is given we 
implicitly take M = M 2 oo, e.g. 



dn 



Hln M 



din M 



n\n m. 



2(H) 



din M vir 
v,r dlnM 200 



(39) 



We also consider the non-linear matter power spectrum 
calculated in the halo model approach (see [13] for a re- 
view). Since all matter is assumed to be within bound 
halos, the matter power spectrum can be decomposed 
into 1-halo and 2-halo terms, 



PmmQt) =I 2 (k)P L (k) + P ih (k) 



(40) 
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10 12 10 13 10 14 10 15 



M 20 „ [M /h] 

FIG. 3: Deviation in the halo mass function at z — of 
sDGP from a dark energy model with the same expansion 
history (QCDM). The points show measurements in the full 
and linearized DGP simulations, the band shows the Sheth- 
Tormen + spherical collapse prediction range between full 
DGP collapse (blue dashed line) and linearized DGP collapse. 




10 ia 1Q13 10 M 10 15 



M 20 „ [M /h] 

FIG. 4: Same as Fig. [3] for the two normal branch DGP + 
dark energy models nDGP-1 (top) and nDGP-2 (bottom) rel- 
ative to ACDM. The linearized DGP simulation results have 
been displaced horizontally for clarity. 



A. Halo Mass Function 



where 

f M 2 
P lh (k) = d\nM vir n lnMvii ^\y(k,M vil )\ 2 , (41) 

J Pm 

/M ■ 
dlnMyir rii n A/ vir ^^y(fc,M vir )&ii n (M vir ) . 

Here, y(k, M) is the Fourier transform of an NFW den- 
sity profile truncated at i? v ir, and normalized so that 
y(k,M) -t 1 as k -> 0. Note that with the ST mass 
function and bias, lim^o I(k) = 1. 



V. RESULTS 

We compare our spherical collapse and halo model 
predictions with the results of N-body simulations pre- 
sented in of the sDGP and nDGP+DE models (see 
§ (Til Tab. |T|. In addition to the full simulations which 
solve the non- linear ip equation simulations using 
the linearized tp equation have been performed through 
Eq. (20! . 

We always compare observables measured in the DGP 
simulations with those of GR simulations with the same 
initial conditions and expansion history. In this way, cos- 
mic variance as well as systematic issues cancel out to a 
large extent. 



Fig. [3] shows the deviation of the halo mass function 
from QCDM measured in the sDGP and linearized sDGP 
simulations, and the spherical collapse predictions. Fig. [4] 
shows the corresponding results for the nDGP+DE mod- 
els. The spherical collapse predictions work well in both 
cases and in particular match the shape of the deviations 
and the relative impact of the Vainshtein mechanism. In 
both cases, they somewhat underestimate the size of the 
deviations at a fixed mass, corresponding roughly to a 
shift in lgM 20 o of - 0.3 - 0.5. 

In the nDGP models, force modifications are stronger 
(see Fig. [2|, leading to larger deviations in the mass func- 
tion from the corresponding GR model with the same ex- 
pansion history. In particular, the abundance of massive 
halos M200 10 14 Mq/H is significantly enhanced. This 
behavior is due to the exponential sensitivity of the mass 
function to v = 5 c /a(M) at the high-mass end, and is 
similar to what was seen in the large- field f(R) models 
in @. 

Furthermore, the density threshold e _1 for the onset of 
the Vainshtein mechanism is higher in the nDGP models 
than the sDGP model, so that the mass function is less 
affected by the Vainshtein mechanism in nDGP. This is 
borne out by both simulations and spherical collapse pre- 
dictions: the relative spread in the predictions between 
the DGP and linearized DGP case shrinks considerably 
when going from large to small r c , i.e. from sDGP to 
nDGP-2 to nDGP-1. 
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FIG. 5: Linear bias bu n of dark matter halos measured in the 
ACDM and QCDM simulations at z = (see g IVBj) . and the 
prediction of the Sheth-Tormen prescription. 



Finally, since the full spherical collapse predictions al- 
ways slightly underestimate the deviations, they can be 
used to place conservative limits on DGP braneworld sce- 
narios from measurements of the halo mass function e.g. 
from massive clusters f35T - l37j . Alternatively, the predic- 
tions can be recalibrated based on simulations by intro- 
ducing a constant shift in lg M200 ~ 0.3 — 0.5. 



B. Halo Bias 

This section presents new results on the clustering of 
halos in the simulations of 0, [26[ . We extract the linear 
halo bias b\[ n (M) from our simulations as described in Q. 
For halos of a given logarithmic mass range in a box of 
size Lbox, we first obtain the halo bias b(k, M) by divid- 
ing the halo- mass cross spectrum Ph m (k) by the matter 
power spectrum for each simulation run. 2 In order to re- 
move trends from the non-linearity of the bias, we then fit 
a linear relation to b(k, M) = b\\ n {M) + a(M) k between 
the minimum k (the fundamental mode of the box) and 
~ 15fc m i n , where b(k,M) is the combined measurement 
from all boxes. The same fitting procedure is applied to 
the run-by-run ratio of &dgp(^ -&0/&GB.(fcj M). 



2 Note that the definition of bias adopted will differ from alter- 
nate choices such as (Phh/^mm) 1 ' 2 or Phh/^hm m the non-linear 
regime where the correlation coefficient between halos and mat- 
ter can differ from unity. 



FIG. 6: Relative deviation of the linear bias Abn n /&ii n in the 
full and linearized sDGP simulations from that of the QCDM 
simulations as a function of halo mass at z = 0. The linearized 
simulation points have been displaced horizontally for clarity. 
The band shows the Sheth-Tormen + spherical collapse pre- 
diction range between full DGP collapse (blue dashed line) 
and linearized DGP collapse. 



We then bootstrap over many realizations of the set of 
simulations, performing the fit for every realization. We 
use the average of the fit parameter bn D (M) as estimate 
of the linear bias, and its spread as an estimate of the 
error. Note that in case of the nDGP simulations, we 
only have 3 runs per box size, so that the error estimate 
itself has a large uncertainty. 

We show the linear bias b\[ n (M) as a function of halo 
mass for the QCDM and ACDM simulations themselves 
in Fig.[5j As the halo mass function deviates significantly 
from a pure power-law, especially at the high-mass end, 
we plot the bias measurements at the position of the mea- 
sured average lg A/200 of the halos. The Sheth-Tormen 
prediction of Eq. (|36[) . using the parameters from Tab. [TT1 
and rescalcd from M v ; r to M200 in each case, matches the 
simulations well for masses up to 10 14 M©//i. At higher 
masses it overpredicts the bias in the simulations, though 
the deviation is a the level of 1 — 2a. Note that halos are 
more biased at a given mass in QCDM than in ACDM, 
due to the reduced growth and smaller power spectrum 
amplitude in this model. 

This trend continues in the DGP simulations: in 
sDGP, gravity is further weakened by the repulsive brane- 
bending mode, so that the linear halo bias at fixed mass 
is increased by ~ 10% compared to QCDM (Fig. [5]). 
There is a hint that halos are slightly higher biased in 
the linearized simulations compared to the full Simula- 




FIG. 7: Same as Fig. [5] for the two normal branch DGP + 
dark energy models nDGP-1 (top) and nDGP-2 (bottom), 
relative to ACDM. 



tions, though both are consistent given the error bars. 
The spherical collapse prediction matches the simulation 
results over the whole mass range. The full DGP collapse 
(blue dashed line) marks the lower edge of the shaded 
band, while the linearized DGP collapse corresponds to 
the upper edge, in accordance with expectations. The 
spread of the predictions is similar in magnitude to the 
tentative differences between linearized and full simula- 
tions. 

For the nDGP models, the opposite trend in the bias is 
seen (Fig. [7]): halos are less biased at a given mass. While 
the full and linearized simulation results are consistent 
in the nDGP models, showing no sign of the Vainshtcin 
mechanism affecting the linear halo bias, they both fol- 
low the spherical collapse prediction very well. This is in 
accordance with the good description of the mass func- 
tion results in § IV Al Also, given the small spread in the 
spherical collapse predictions, we do not expect to see 
any differences between full and linearized simulations 
for the nDGP models. 



C. Non-linear matter power spectrum 

We can now assemble the ingredients of the halo model: 
the mass function, bias, and profile of halos, to predict 
the non- linear matter power spectrum following §HV] We 
assume that the inner parts of the density profiles of halos 
are not affected by DGP, as indicated by simulations (2(| . 
More precisely, one can assume that for a halo of given 
mass M200, the scale radius r s (Eq. ([57)) ) is the same 



FIG. 8: Deviation in the matter power of the sDGP model 
from QCDM at z — 0. The points show measurements in the 
full and linearized DGP simulations, while the band shows 
the halo model prediction based on spherical collapse and the 
Sheth-Tormen prescription (between DGP [blue dashed line] 
and linearized DGP collapse). The long-dashed line shows 
the renormalized perturbation theory prediction from [30I ]. 



in DGP as in GR. Adopting Eq. (|38[) for the concentra- 
tion c v ir = Rvir/fs m GR, we have for the concentration 
relation in DGP: 



Cvir, DGP 



(M) 



Ra(M) 



A vir (DGP) 



Ra(M) 



:(M vir , GR ), (42) 



A vir (GR) 



where R A (M) = (3M/(4nAp))^ 3 , and M vil , GR is the 
virial mass in GR that corresponds to a virial mass of 
M in DGP. We find that for all DGP models considered 
here, the concentration c v i r ,DGP defined by Eq. (|42[) is 
within 3% of the standard relation c v i r (M V i ri DGp)i which 
has a negligible effect on the power spectrum on the scales 
probed by the simulations. Hence, we leave the concen- 
tration relation Eq. (|38|l unchanged in our power spec- 
trum predictions. 

Fig. [8] shows that the sDGP simulation results are 
matched very well: the DGP collapse prediction (blue 
dashed line) is close the to full simulation results, while 
the linearized DGP collapse is close to the linearized sim- 
ulations. The renormalized perturbation theory predic- 
tion of [3p| which uses a completely different approach 
to take into account the non-linear interactions of the 
brane-bending mode is also quite close to our DGP col- 
lapse calculation. 

The match to the power spectrum in the nDGP simu- 
lations (Fig. [5]) is somewhat worse, showing discrepan- 
cies of - 30% for nDGP-1 and - 10% for nDGP-2. 
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FIG. 9: Same as Fig. [5] for the two normal branch DGP 
+ dark energy models nDGP-1 (top) and nDGP-2 (bottom) 
relative to ACDM. The long-dashed lines show the predictions 
of halofit using the linear DGP power spectrum. 



In particular, discrepancies can be seen in the quasi- 
linear to non-linear transition, k ~ 0.1 — 1/i/Mpc for 
nDGP-1. The simplified l-halo/2-halo split in the halo 
model breaks down on these scales. We also show the 
predicted non-linear power spectra using halofit [38| 
in combination with the linear DGP power spectra at 
z = 0. Note that for k > 0.01 h/Mpc, the linear DGP 
power spectrum is identical to that of a ACDM model 
with higher linear normalization. While the match to 
the simulations is better than our spherical collapse pre- 
dictions at k < 0.1/i/Mpc, the deviations grow towards 
smaller scales so that halofit is a worse fit to the simula- 
tions than the spherical collapse model for k >1 h/Mpc. 
Clearly, it is not trivial to model the simulation results 
for AP(k)/P(k) to better than - 10%. 



While the overall magnitude of the power spectrum 
enhancement in nDGP is not matched, the shape of 
AP(k)/ P(k) is matched quite well. Furthermore, the ef- 
fect of the Vainshtein mechanism on the power spectrum 
is predicted accurately, as can be seen by comparing the 
spread in the spherical collapse predictions to the differ- 
ence between full and linearized simulation results. These 
findings indicate that it should be possible to rescale lin- 
earized DGP simulation results with the relative Vain- 
shtein suppression calculated from spherical collapse. 



VI. CONCLUSIONS 

By studying the collapse of a spherical perturbation 
under DGP braneworld gravity, we have shown how sim- 
ulation results on the mass function, halo bias and power 
spectrum can be understood semi-analytically. In DGP 
gravity, force modifications are carried by the scalar 
branc-bending mode. The global properties of the re- 
sponse of the brane-bending mode to matter control how 
the Vainshtein mechanism modifies force and energy con- 
ditions during collapse. These conditions are important 
for the calculation of virial equilibrium (see the Appendix 
for detailed discussions of these results). 

In particular, the presence of evolving modifications to 
the gravitational force either through the brane bending 
mode or through the background expansion violate con- 
servation of Newtonian total energy for traditional defini- 
tions of the potential energy contribution. This violation 
applies to smooth dark energy models with w ^ -1 as 
well. Wc introduce a new, general technique for defin- 
ing the virial radius which does not rely on strict energy 
conservation. 

Under the halo model, these spherical collapse predic- 
tions give rise to predictions for the mass function, halo 
bias and power spectrum. We have shown that these 
predictions are in good qualitative agreement with DGP 
N-body simulations on both the self-accelerating and nor- 
mal branch. In particular, the use of spherical collapse 
for the mass function always provides slightly conserva- 
tive limits on mass function deviations when compared 
with the simulations. Hence the semi-analytic techniques 
introduced here can be used as a practical tool for extend- 
ing simulation results for the purpose of studying param- 
eter constraints on braneworld models from observations 
of the mass function. 

While the absolute power spectrum agreement is not 
quite as good, these techniques can still be useful in com- 
bination with linearized DGP simulations. Since our 
spherical collapse predictions appear to capture accu- 
rately the impact of the Vainshtein mechanism on the 
power spectrum, they could provide an effective way of 
taking non-linear interactions into account in results ob- 
tained from linearized DGP simulations. Such simula- 
tions are very easy to implement and an order of mag- 
nitude cheaper computationally than the full DGP sim- 
ulations. Likewise they are more readily extendable to 
higher resolution and can be used to cover a wider range 
in parameter space. 
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Appendix A: Top-Hat Overdensity in DGP 

In this Appendix we present in detail the various as- 
pects of top-hat perturbations in DGP used in the main 
text. We begin by reviewing the techniques introduced 
in Ref. [f| for the brane bending field but pay special at- 
tention to the matching of the various solutions as well 
as derive a closed form expression for the global field pro- 
file. This global profile is used to study virial equilibrium 
( QA 3[) and potential energy ( £1A 4[) . We discuss subtleties 
due to the relationship between the two and violations of 
Newtonian energy conservation in QA 51 and their impact 
on defining the virial overdensity in §A 61 Finally we ex- 
amine the impact of density compensation in the exterior 
of the top-hat in $A 71 in light of the lack of a Birkhoff 
theorem in DGP. 



1. Spherical Symmetry 

We begin by assuming that the density perturbation 
5p(r) is spherically symmetric but otherwise arbitrary. 
Given the induced spherical symmetry in the field solu- 
tion, the field equation © reduces to 



1 d 2 d<p ^t c 
r 2 dr dr 3/3 

8ttG 



4 d 2 ip dip 
r dr 2 dr 



1 dip 
r dr 



3/* 



■Sp. 



(Al) 



Integration of the field equation (|A1|) over r 2 dr then 
yields 



,dip 



2 2 ( dt f^2 



-Gm{r), 



dr 3fi c y dr' 3/3 
with the enclosed mass perturbation defined as 

m(r) = 47r / r' 2 Sp(r) dr' . 
Jo 



(A2) 



(A3) 



If 6p(r) = for r > R, the enclosed mass fluctuation 
m(r > R) is the total mass fluctuation SM. 

Note that the ViVjtp term in Eq. © is critical in ob- 
taining this solution since it causes cancelation of the 



integral terms when integrating by parts, leaving only 
the boundary terms for the non-linear piece. The field 
solution interior to r has no direct effect on the solution 
at r. Like Newtonian dynamics, only the enclosed mass 
not the enclosed field matters. This property is crucial 
for maintaining the linearity of the field solution at large 
r in the presence of strong non-linearity at small r. More 
generally, linearity is a consequence of Eq. Q satisfying 
a non-linear Gauss's law and does not require spherical 
symmetry (see e.g. j39j). 

Since m — > SM at distances beyond which there are no 
density fluctuations, and assuming that <p(oo) = 0, we 
can immediately see that the far exterior solution must 
be to leading order tp oc 1/r. Given an increasingly small 
non-linear term, this requires 



lim ip = — 



2 GSM 

3/3 r 



(A4) 



For the small r solution, the key simplification is that 
the force modification dtp/dr is now an analytic function 
of the enclosed mass 



dp 3/3r*(r) 



dr Ar 2 

where the Vainshtein radius of the enclosed mass is 

, . fl6Gm(r)r 2 ^ 1/3 
r * {r) ={ 9/32 



(A5) 



(A6) 



and 



g{x) = x[\J 1 + x 3 — 1] 



(A7) 



Note that in general r* is a function of r and reflects the 
enclosed mass, and r 3 /r 3 the average density, not the 
local density. The latter would be implied by setting 



(ViV^) 2 = c (VV) 2 



(A8) 



with c =const. in the original field equation. However, 
doing so would violate the r — > oo limit. In particular, 
the far field limit in this approximation would reveal the 
presence of a Vainshtein screened mass instead of the true 
mass perturbation 5M . Small scale non-linearity in the 
density field would then no longer average to give the 
required linear perturbations on large scales [40| . Such 
an approximation or any that relates the field solution 
to the local density should not be used in a cosmological 
context (cf. [Mill)- 



2. Field Profile 

We now assume a top-hat spherical perturbation 
of radius R with a constant density enhancement Sp 
[Eq. dTTJ)] . As usual we neglect any compensating under- 
density swept out by the prior evolution of the top-hat. 
Since the force modifications within the top-hat only de- 
pend on the enclosed mass m(r) and not the exterior, 
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the compensation does not impact the dynamics. It can 
however influence the field profile r > R and we return to 
this point in 71 For the pure top-hat, m(r > R) = SM 
and we can solve for the profile tp(r) in closed form. 

First let us consider the exterior solution at r > R. In 
the exterior m=const. and there is a single Vainshtein 
radius R* = r*(R). Defining a new variable x = r/R*, 
we can write Eq. (IA5|) as 



where 



A = 



dp 
dx 



3/3 (R 



Ag(x), 



(A9) 



4 GSM 



4 \r c J 3/3 R, ' 
and we can obtain the full exterior field solution as 

dip 



(A10) 



ip(r) 



r/R, 



dx 



-dx 



(All) 



Ax 2 



{2^1-1/2,-2/3,1/3,-1/a; 3 ]-!} 



The solution of course recovers Eq. (|A4[) in the limit of 
x > 1 



2 GSM 



lim p = : 

*»i 2x 3/3 r 



(A12) 



Note that the constant A is — 2ip* where tp* is the lin- 
earized field profile evaluated at the Vainshtein radius. 
In the opposite limit x <C 1 



lim p = -A(C - 2x 1/2 ). 



(A13) 



The constant piece, 



Co 



A~ 



ax— 

dx 



r[i/3]r[i/6] 



4^ 



therefore dominates, and 



lim w : 



4C GSM 
3/3 R* 



2.103, 



= const. 



(A14) 



(A15) 



However field gradients are determined solely by the x 1 ^ 2 
term. Since particle dynamics depend on forces and 
hence field gradients, the existence of a constant term 
in addition to the x 1 ^ 2 term in the field profile makes an 
important difference when comparing energy conditions 
like conservation laws and dynamical considerations like 
virial equilibrium, as we shall see. This distinction is of- 
ten neglected in the literature (e.g. [H, HI])- Fig- E3 
shows \ip(SM; R) \ in units of the linearized (x ^> 1) value 
as a function of R/R*. Note the strong suppression of 
the surface potential and its saturation for R <C R* ■ 
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FIG. 10: Surface field profile \ip(SM, R)\ (solid curve, see 
Eq. dAT8j) in units of (2/3/3) GSM/R*, the linearized value 
at i?*, as function of xr — R/R t . Parameters are for the 
sDGP model (a = 1). The dotted line shows the linearized 
solution, xT, 1 in these units. In the Vainshtein limit xr <C 1, 



ip(8M,R) approaches a constant, 
the binding energy U v (S IA 4[> . 



p{8M, R) also determines 



The field interior to the top-hat can be similarly solved. 
Since r* (r) cx r for a constant density profile, Eq. (|A5j) 
implies dip/dr oc r and hence 



p(r) - p(0) = Br 2 



We see from substitution into Eq. (|Al[) that 



B 



A 



2RR 



(A16) 



(A17) 



solves the field equation. Here xr = x(R) = R/R*. 
Note that d<p/dr is automatically matched at r = R to 
the exterior solution Eq. (| A9|) . The remaining condition 
is that the interior field solution at r = i?, 



<p(0) 



(A18) 



match the exterior solution from Eq. (| Al 1|1 . Fig. [TT] 
shows the resulting ip profile for different values of xr, in 
comparison with the Newtonian potential "J at. Again it 
is interesting to examine the x R 3> 1 and x R <^ 1 limits 
of Eq. (|A18[) . In the former case we regain the linearized 
(Newtonian) expectation that the central value is 3/2 of 
the surface value 



lim ip(0) 

XR>1 



GSM 
PR 



(A19) 
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FIG. 11: Field profile f(r) for a top-hat mass of radius R 
in units of the Newtonian surface potential G8M/R, for two 
values of the Vainshtein radius R*/R = 10, 100i?. if was 
scaled by 3/3/2 so that the linearized field solution agrees with 
the Newtonian potential shown as the long-dashed line. 



In the opposite, "Vainshtein" limit, the central field like 
the surface field is independent of R to leading order 



ACqGSM 

lim V «p(J2 . 
zi?«i dp R* 



(A20) 



The central field value, which determines the potential 
energy associated with the fluctuation is therefore sup- 
pressed by (4Cq/3)R/R* from the linearized expecta- 
tion. Moreover the change in the field interior or near 
the object, important for forces, is further suppressed by 
(i?/i?*) 1/2 . We discuss the consequences of these features 
of the field profile for virial equilibrium and the potential 
energy in the following section. 

The force suppression can be recast in terms of an ef- 
fective Newton constant in the Poisson equation for the 
ip field. Note that in the interior solution, the two pieces 
of the nonlinear term in Eq. ^ combine to form 



(VV) 2 - (V 4 V^) 2 = x(VW ■ (A21) 



Since the field equation is then algebraic in V 2 y>, one can 
solve for V 2 ip to obtain: 

VV = SttAG BG p(R/R*)Sp, (A22) 
2 

AG DGP {x) = —g( x )x 2 G 
2 

= 3^ 



[Vl + x~ 3 -l]x 3 G. (A23) 



Hence, in combination with the Newtonian piece, the in- 
terior solution can be phrased as possessing an effective 
Gdgp = G + AGdgp modification to gravity. Further- 
more x~ 3 cx S and so Gdgp is a function of the local 
density inside the top-hat. This relation is specific to 
the interior of a top hat and is not simply a spherically 
symmetric approximation as the exterior solution shows. 
More generally, we can see by taking the derivative of 
Eq. (lATfll that 



2Gm(r) d 2 



VV = 87rAG DGP (r/r.) 6p(r) + 

r~ ar 

(A24) 

Note that AGdgp, m(r) and r*(r) are not local functions 
of the density field but involve the full interior profile. We 
shall return to this point in £|A 71 



3. Virial Equilibrium 

The virial theorem arises from integrating over space 
the first moment of the Boltzmann equation, i.e. from the 
equation of momentum conservation (see e.g. [4l|, §4.3). 
Despite its usual association with potential energy, the 
virial theorem is inherently a force balance equation and 
is the collisionless analogue of hydrostatic equilibrium. 
Thus the virial condition is immune to ambiguities in 
the definition of potential energy that we shall discuss in 

El 

The virial theorem reads 



W 



d 3 x p m (x) x- V* = -2T, 



(A25) 



where W receives contribution from the Newtonian grav- 
ity of the overdensity, the effective background term, and 
the brane-bending mode tp. For a spherically symmetric 
top-hat each contribution to \& yields 



IT"; 



-3M 



r 2 dr d^i 
~W T ~dV 



(A26) 



where VE^ stands for either VP jv , ^eff, o r = Thus, 
any constant offsets in "1^ do not contribute to W. In 
particular the constant term in ip in the Vainshtein limit 
of Eq. (|A20p and its implied potential energy does not 
enter into the virial condition. 

In the case of the Newtonian contribution, Wn defined 
in this way is given by: 



W N 



3 GMSM 
5 R ' 



(A27) 



For the effective contribution of the background ^ ff cx r 2 
and 



W cS 



4ttG 



[(l + 3w eS ) PcS + p}MR 2 . (A28) 



Note that in our convention, we have included the p term 
in Weft rather than Wn- Adding the two contributions 
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yields the familiar result [21] : 

3 GM 2 AttG 
5 _ R 5~~ 



W N +W ( 



cir 



{l+3w cS ) PeS MR 2 . (A29) 



Since the trace of the potential energy tensor is defined 
via forces, the ip contribution can be described in terms of 
AGdgp- First let us examine the exterior region. Using 
Eq. (1X51). we obtain 



dip 
dr 



4 GSM 

3/3 r 2 



x 2 g(x), 



which using Eq. (|15[) becomes 

d^ = AG DG p(r/R*)SM 
dr r 2 



(A30) 



(A31) 



In the linear regime r ^> R*, AGdgp then reduces 
to (3/3) _1 G. Note that in the cosmological context, 
AGdgp also has a slow time dependence through /3(a) 
[see Eq. ([TO])]. 

Using that in the interior dip/dr oc r, we have 



3M 



3M 



10 rflnr 



(A32) 



Given that dip/dr at r = R is determined by the exterior 
solution, we obtain 



3 AGdgp A/fl A/ 
5 # ' 



(A33) 



where AGdgp = Gdgp — G is the effective gravitational 
constant for the force modification [Eq. ([I5])]. Adding all 
three contributions, we obtain Eq. (|31[) . 

Specifically in the collapse calculation, we evaluate the 
virial condition Eq. (|A25I) in terms of the kinetic and 
potential energies per unit mass, written in terms of our 
collapse variable y(lna). Defining 



E = ^M{H Rif 



(A34) 



where Ri is the initial radius of the perturbation, we 
obtain 



H 2 

— j(y' + a/a t ) 2 , 



— - = -a m a- 3 (5 + l)(y + a/ ai y 



T 
Eo 

W_N 
Eq 



; AGugp{R/R* 
G 



■S (y + a/a. t ) 2 , 



= -[1 



3K; e ff(a)] Peff ° (y + a/ai) 2 . 



Pcr.O 



(A35) 



We then define the virial radius as the radius during the 
collapse at which the virial condition is satisfied. We shall 
examine approximate techniques for finding this scale 
through energy conservation in 61 As we shall see, 
with conventional definitions of potential energy, the to- 
tal energy is not strictly conserved, especially during the 
initial stages of collapse. 



4. Potential Energy Definition 

Let us now consider the potential or binding energy 
of the top-hat mass. For the Newtonian contribution, 
the potential energy is well defined. By virtue of the 
Birkhoff theorem, we can view total mass inside the top- 
hat as a Newtonian system in a flat background. Wc 
shall see that neither the potential energy contributed 
by the brane bending mode <p nor the effective forces of 
the background expansion are unambiguous to define as 
both depend on the exterior cosmological context. We 
therefore follow the convention in the literature in defin- 
ing them by analogy to the Newtonian contribution. 

The Newtonian calculation of potential energy pro- 
ceeds by replacing the exterior of the top-hat with a 
flat background, the metric analogue to the Newton iron 
sphere theorem (see [42j §4). Removing mass shell by 
mass shell from the outside in, we obtain 



U = 



^tot(m(r);r) 4irp m r 2 dr 



(A36) 



where ^t t(m; r) denotes the solution for the total grav- 
itational potential for a top-hat with radius r and mass 
perturbation m = (r/R) 3 SM evaluated at r. This is not 
to be confused with the total potential at a radius r in- 
terior to the whole top-hat. 

Even for the matter contribution there arises an ambi- 
guity in that the contribution of the background matter 
density across the top-hat is defined only up to a con- 
stant through ^ e ff of Eq. P^]) . In other words, the 
iron sphere theorem applies directly to forces not poten- 
tials and constant offsets do not have any impact on the 
dynamics. Therefore 



N 



off 



*0- 



(A37) 



This constant does not simply introduce a trivial shift 
in U since even though it is constant across the top-hat, 
it is not necessarily constant as we strip away mass shells. 
In fact, it is conventional to choose this value to corre- 
spond to the result from Newtonian mechanics 



3 GM 



= -2irGpR z = -- 



2 R 



(A38) 



such that 



GA1 '2 .t ( i 

*w + *cff + *o = ^- + :L -^{l + 3w cS ) PoS R 2 . (A39) 

The —GM/R piece then corresponds to the Newtonian 
mechanics result for the potential given the total mass 
inside the top-hat. Using the Birkhoff theorem, this is a 
valid interpretation of the cosmological case as well. 

Note that the p e ff piece cannot properly be consid- 
ered a binding energy since its contribution cannot be 
considered without reference to the cosmological back- 
ground. Even for quintessence models where p e s rep- 
resents a real energy density, this contribution is sup- 
posed to be smooth within its horizon sized Jeans scale 
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regardless of the top-hat collapse and so excising the top- 
hat and placing the mass in a flat background does not 
strictly make sense. 

Nevertheless, under this convention the binding energy 
from these three components becomes 



u N +u eS +u 



3 GM 2 2ttG 
5 R 



{l+'3w cS )MR 2 . (A40) 



By analogy let us compute the potential energy contri- 
bution from 



1 

2 



tp(m(r);r) Airp m r 2 dr. 



(A41) 



Again, tp(m; r) denotes the exterior solution of tp from 
Eq. (|A11[) for a mass SM — > m and radius R — > r. Note 
that by making this assumption we are implicitly invok- 
ing the Birkhoff theorem where it does not in fact strictly 
apply. Specifically, as mass shells are stripped away, we 
ignore the impact of the underdensity left behind outside 
of the body. As we shall see in ^A 71 this underdensity 
actually changes the interior profile. Nonetheless, since 
the ambiguity mainly affects the initial stages of collapse 
when SM <C M, it is useful to simply define Eq. (|A41[) 
as the binding energy associated with tp for energy book- 
keeping purposes. 

Re-expressing in terms of the dimensionless radius y = 
r/R (not to be confused with y defined in § IIIII which 
we shall not use hereafter), we have m = y 3 M, so that 
x r = r/r*(m) = R/R* = xr is invariant, while A(m) = 
y 2 A(5M). Therefore, 



tp(m;r) = y 2 tp(SM;R), 



(A42) 



where tp{SM; R) is the solution for the full mass eval- 
uated at r = R [Eq. (|A11|) at x = xr\. This is the 
same scaling that one obtains for a Newtonian potential, 
fy(SM;R) = GSM/R and likewise the potential energy 
follows the same scaling 



U v = -Mtp(SM;R) / y*dy 
z Jo 

= ^Mtp(SM;R). 



(A43) 



The important difference for U v is that ip(5M; R), the 
field profile at R, behaves very differently in the linearized 
and Vainshtein limits 



U v = 



1 3 GMSM 
3/3 5 R ' 

2Ca GMSM 



5/3 



xr > 1, 
Xr < 1. 



(A44) 



Hence, unlike the Newtonian binding energy, \U V \ has 
a maximum value which is reached asymptotically as 
xr = R/R* decreases. Thus, while the Newtonian bind- 
ing energy can become an arbitrarily large fraction of 
the rest mass energy SM, the energy in tp is limited to a 
fraction of GSM//3R* ~ (GSM/r c ) 2 ^ 3 due to the brane- 
bending mode interactions regardless of the value of R. 



5. Potential Energy Usage 

Defining a Newtonian based potential energy even 
though the collapse does not require a Newtonian inter- 
pretation is useful for two interrelated reasons. Firstly it 
serves as a bookkeeping device if the total energy is con- 
served during the collapse. Secondly, it can be used to 
evaluate the virial condition if it can be simply related to 
the trace of the potential energy tensor W. We examine 
to what extent these two expectations are satisfied given 
the field tp and the effect of the background expansion. 

Let us define the total energy of the perturbation dur- 
ing collapse as 



E = T+U. 



(A45) 



Taking the time derivative of Eq. (|A45|) and using the 
equation of motion of R(t) [Eq. (|18p] rewritten using the 
top-hat profile as: 

n = _9i^_^ [il + 3Wes)PeS+m (A46) 



we obtain: 



R 2 



dE 3 •» dU ■ dU 
H - 5 MRR+ 0R R+ lH 

dU 

dt 1 



(A47) 



Here, the partial derivative dU/dt receives contributions 
from evolving quantities in the total potential energy, the 
sum of Eq. (|A40j) and (|A43|I . Since M is conserved, there 
is no violation of energy conservation for a pure matter 
system with p e g = and tp = 0. This is a consequence 
of adding the \l/o offset term to reproduce Newtonian 
mechanics in the matter terms including the background. 

On the other hand both the p e g and tp contributions 
have explicit time evolution across a Hubble time. Note 
that violation of energy conservation due to evolution of 
Peg also applies to dark energy models where w e g 7^—1. 
In the tp term given by Eq. (|A43[) . f3 evolves with the 
expansion, and as long as S is not much greater than 1, 
SM also evolves during collapse. More generally these ef- 
fects occur whenever the modification to the background 
(due to modified gravity or the presence of dark energy) 
does not match that of the perturbations. We return to 
the impact of energy non-conservation in § IA 61 Non- 
conservation of the Newtonian total energy defined in 
Eq. (|A45[) docs not mean a violation in covariant conser- 
vation of energy and momentum. 

In order to use the potential energy to assist the evalua- 
tion of the virial theorem, we must relate U to the trace of 
the potential energy tensor W. It is well known that for a 
potential satisfying *5?(r) cx r a and for which the interior 
solution is * int (r < R) = * cxt (m(r); r), W t = -aU l . 
Hence, for potentials satisfying this condition, the po- 
tential energy determined by the potential itself is of the 
same order as the trace of the potential energy tensor 
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defined by the forces. This holds for the Newtonian con- 
tribution, where we have Un = Wat, and the effective 
background contribution, which satisfies U e s = — W e g/2. 
For the brane-bending mode, the potential is no longer 
a pure power law and the distinction between U and W 
leads to interesting consequences in the Vainshtein limit. 

Note that in this xr <C 1 limit, the potential energy U v 
is dominated by the constant term in Eq. (|A13[) . while 
the contribution from the x 1 ^ 2 part of the profile which 
determines forces is much smaller. Correspondingly, the 
assumption that the trace of the potential energy tensor 
W is of order the potential energy U is not valid for the 
if contribution, as the change in the potential across the 
body is much smaller than the potential depth itself. The 
relationship between W v and U v follows from Eq. (|A32j) 
and ([A43]) : 



U 



r=R 



d In Up 



(A48) 



and W v = U v as usual 



In the R 3> -ff* limit ip oc r 
but in the R <C R* limit, ip ~ const, and the trace 
of the potential tensor is highly suppressed compared to 
the potential energy. Fig. [T^] shows U v and W v as a 
function of the overdensity 8 = 5p/~p of the perturbation. 
Note that if we were to interpret the Vainshtein effect 
as simply a modification of G we would infer the wrong 
energy condition at virialization. We discuss this issue in 
the next section. 
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FIG. 12: Scaled brane-bending mode binding energy 3/3 x 
U v in units of the Newtonian binding energy Un = 
—3GM5A1/5R (black solid) as a function of the overdensity 5 
for a spherical top-hat mass at z = in the sDGP cosmology. 
The blue dashed line shows the trace of the potential energy 
tensor (3/3) x W v used in the virial theorem, Eq. (|A32[) . again 
with respect to the Newtonian value Wn = Un- 



6. Misestimating Virial Overdensity 

Most commonly, the virial condition Eq. ([29]) is evalu- 
ated using energy conservation. At virialization, 

T(i? vir ) = -\w{Ry lt ) (A49) 

and the total energy E = T(R vir ) + U(R vir ) = U(R vil -) - 
W(R v i r )/2. Since at turn-around (R = i? ta ) the kinetic 
energy vanishes, we have E = [/(i? ta ). Assuming energy 
conservation we obtain 

C/(i? ta ) - U(R vh .) + \w{R vn ) = 0. (A50) 

By further assuming a relationship between the potential 
energy U and the trace of the potential tensor W we can 
solve for i? v ir/ Rt&- There are therefore two ways by which 
this association can go wrong: if E is not conserved and 
if an incorrect relation between U and W is employed. 

Let us begin by examining the first issue. Energy is 
not strictly conserved in any model where either w e fi ^ 
— 1 (including quintessence), modifications to gravity are 
time variable, or force modifications are only generated 
by the perturbed mass SM. 

Evaluating Eq. (|A45|) during collapse, we found that 
in the effective dark energy model QCDM (which has 
the same H(z) as sDGP), energy conservation is vio- 
lated by ~ 3% from turn-around to collapse. While the 



violation of energy conservation in dark energy models 
thus seems to be minor, it does influence the virial over- 
density due to the sensitivity of A v j r to i? V j r and a V j r . 
Fig. [13] shows A v j r as function of a constant dark energy 
equation of state w e g and fl e s = 0.741, determined by 
evaluating the virial condition during collapse (our ap- 
proach) and using energy conservation (e.g., [2l|). While 
both approaches agree for w c s — — 1 as expected, there 
arc clear differences as soon as w e s ^ — 1. Note that 
these differences are of the same order as the difference 
A vir (w cS ) - A vil (w cS = -1). 

Our approach does not rely on exact energy conser- 
vation and one might thus expect the A v j r obtained in 
this way would lead to a better match to observablcs. 
Since quantities such as the mass function are typically 
simply fit to simulations with a given definition of over- 
density, this is in part an issue of semantics. However 
use of a more physically motivated scaling might lead 
to a more universal form for the mass function or one 
that scales more simply with parameters of the theory. 
It would however be necessary to compare with N-body 
simulations of w c s =/= — 1 dark energy models and that 
is beyond the scope of this work. Here we simply note 
that the dependence on w e s of A V j r determined using our 
approach is smaller than that of the usual definition and 
so it would predict a more universal scaling with a fixed 
overdensity than the standard approach. 
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FIG. 13: Misestimation of the virial overdensity A v i r in 
quintessence dark energy models. A v ir as defined in § IIII CI 
for collapse at oo = 1, as a function of the dark energy equa- 
tion of state parameter u> e ff =const. and fi e ff = 0.741 i.e. 
for standard gravitational forces (top panel). The solid line 
shows Avir determined by evaluating the virial condition dur- 
ing collapse, the approach adopted here, while the dashed line 
shows the usual calculation using standard energy conserva- 
tion as described in S IA 61 Both calculations agree for w = — 1, 
i.e. p c fj =const. The lower panel shows the relative deviation 
between the two, (A v " - A v ir)/A V ir. 



The deviations from E = const, are much larger in 
the DGP modified force case, where evolving forces and 
6M lead to stronger evolution of E. Here, differences in 
the total energy between turn-around and collapse are 
10 — 15%. Thus, the choice of procedure for determining 
A v i r becomes even more important. Again, we found 
that our approach (greatly) reduces the dependence of 
A V i r on the evolution of the modified forces. 

Finally, we consider the effect of assuming an incorrect 
relation between U and W in Eq. (|A50|) . For example, 
it is tempting to just set U v = W v = AGdgp/G x Wn, 
i.e., ignoring the R- independent term in U v . However, 
since AGdgp — > in the late stages of collapse (S — > oo), 
U v would erroneously be set to in this approximation. 
For /3 of order unity, this leads to apparent violations 
of energy conservation at the level of 30%. When using 
energy conservation in the presence of the tp field, it is 
thus crucial to take into account the differences between 
Uu, and W w . 



FIG. 14; Field gradient dip/dr for a compensated top-hat 
profile (solid) in units of 2/(3/3)GSM/R 2 as a function of y = 
r/R. y e = R e /R is set to 100 (corresponding to S = 10 6 ), 
and y, = R* / R — 67, corresponding to t = 0.3 valid for 
the sDGP model. The dashed line shows dip/dr for the same 
top-hat mass but with uncompensated profile. Near R, force 
modifications from <p are not affected by compensation. 

7. Compensated Top-hat Profile 

Finally we study one example of a density profile be- 
yond the pure top-hat, the compensated top-hat profile. 
In the cosmological context, a collapsing top-hat pertur- 
bation sweeps out "empty space" and in fact has the 
following density profile; 

f pS, r< R, 
p(r)-p = < -p, R<r<R e , (A51) 
[ 0, r>R e , 

where R e = Ria/ai is the physical radius today corre- 
sponding to the radius of the perturbation at an early 
time a% when S -^i 1. We continue to call the total 
mass and mass perturbation enclosed at R as M and 
SM = m(R) respectively. 

In terms of the scaled radial coordinate y = r/R, the 
full description for the enclosed mass perturbation be- 
comes 

( SMy 3 , y<l, 

m(r) = <^ M[l - {y/Vef], Ky<y e , (A52) 

Importantly, forces at a given radius only depend on the 
enclosed mass m(r) through Eq. (|A5[) . Given that for 
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r < R the enclosed mass of the compensated and un- 
compensated top-hat are the same, compensation has 
no impact on the interior dynamics of collapse. Like- 
wise, from the definition of the potential energy tensor 
Eq. (|A26|) . we see then that W v is unchanged from that 
of the pure top-hat profile and hence the virial condition 
is unmodified. 

Naively, one might assume that as long as 8 1, the 
compensation has little effect on the binding energy or 
gravitational potential ^ in the interior as well but we 
shall see that this is not necessarily so for the brane- 
bending contribution due to the Vainshtein suppression. 

In the exterior (y > 1) forces from tp are modified by 
the compensation as 



dtp 
dy 



R 



d(p 
dr 



(A53) 
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In the limit y e 3> y > 1, the forces are the same as in the 
uncompensated profile (see Fig. [14]). Since y^ = 1 + S, 
forces are unchanged near the body, y ~ 1, as long as S 3> 
1. Even for 5 < 1, this modification does not introduce 
any physical effect on the collapse since there is no mass 
in the exterior region 1 < y < y e that could be moved by 
the modified forces. Hence an initial top-hat profile will 
remain a top-hat during collapse. 

Now let us look at the effect of compensation near the 
Vainshtein scale of the mass, = R r /R. Note that 



, 3 e i + s 



(A54) 



where e~ x is the density threshold beyond which the 
Vainshtein mechanism operates as defined in Eq. ([2"T|). 
Unless this density threshold also satisfies e _1 1, com- 
pensation effects will change how the profile saturates, 
since y* will be comparable to y e . Correspondingly in 
Eq. (|A53[) . both S 1 and e <C 1 are necessary for 
dip I dy to recover the uncompensated result at the Vain- 
shtein scale y*. Given that the top-hat ip profile within 
i?* is controlled by its value at i?* , we expect the ip pro- 
file itself to be modified near the body by the density 
compensation unless e< 1. 

More specifically, let us consider the linearized (y» <C 
y) and Vainshtein (y* ^> y) limits as before. First note 
that even the Newtonian force d^pf/dr for the compen- 
sated top-hat at y > 1 is modified as 
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In the linearized limit, Eq. (|A53|) reduces to 
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FIG. 15: The exterior ip profile in units of 2/(3/3)GSM/R 
as function of the scaled radius y for a compensated top-hat 
profile (solid). y e — 100 and = 67 as in Fig. 1141 The field 
goes to at y = y e (dotted vertical line). The dashed line 
shows the profile for an uncompensated top-hat with the same 
mass and radius. Note that the tp field for the two profiles 
differs at all radii, reflecting the fact that tp(r) is not simply 
determined by the enclosed mass at r. 



as expected. In the Vainshtein limit of 3> y(> 1), and 
in which case y <C y e as well, the force contributions take 
the form 



dip 
Inn — — = 



2 G8M\ 2 _ 



y<<y, dy \3f3 R ') ^y V * \y*-l 

(A57) 

We thus recover the leading r 1 / 2 behavior of the force 
in this limit. Fig. 1141 shows dtp/dy as function of y for 
different values of y* = R* /R. Note that unlike the New- 
tonian case, the force given by dtp/dy differs from that 
of an equivalent top-hat with radius r and enclosed mass 
fluctuation m(r) for any r > R. This behavior cannot be 
described by a simple Gdgp(^) parametrization. 

Next as with the pure top-hat, we can solve for the 
whole field profile given the boundary condition tp(R e ) = 
by integrating 
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1/2 



tp(r > R) 
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dp 
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(A58) 



We then obtain 

2 G5M\ y| 



<p(y) = - 



where 
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V --l + F{l,e)-F(lL,e) 
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F(x,e) =4^2^! [-1/2, 1/6; 7/6; ^(l-e" 1 )]. (A59) 
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Note that this is a slightly different hypcrgcomctric func- 
tion than that in the pure top-hat field solution. Fig. [T5l 
shows ip(y) vs. y for a fixed overdensity S = 10 6 for two 
different values of y* (or, equivalently, e). Again, we ob- 
tain the expected scaling in the limiting cases. For the 
linear limit y 3> y*, we have 



lim p(y) 



2 GSM 
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JV- 
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25 y<5 + 2 <5 

(A60) 



proportional to the Newtonian potential for a com- 
pensated top-hat with the same boundary condition, 
^n(Rc) = 0. As expected, for 6 ^s> 1 the profile matches 
the uncompensated case until y approaches y e . 

In the Vainshtein limit y* « e 1 ^ 3 y e ^ 1 and the profile 
becomes 



lim ip(y)-- 
y<y* 



4 GcMA 
3/3 i?» ) 



C £ - 2^/yZ (A61) 



where 



C, 



2e 2 / 3 ' 



(A62) 



In the limit e — > 0, C e = Co, and the profile returns to the 
uncompensated form of Eq. (|A20|> . In Fig. [if)] we show 
C £ /Co the reduction in the surface potential tp(R) due to 
the compensation in the Vainshtein limit. 

One can also define an alternate definition of the bind- 
ing energy of the perturbation SM. Suppose that we 
again define the binding energy U v as in Eq. (|A41[) by 
removing shell by shell of the mass to R e . However, in 
this case we properly account for the impact of the exte- 
rior. As each shell is removed to i? e , it fills in the mass 
deficit so that R e decreases in such a way as to keep 
y e = R e /R constant. 

Thus, the only part of Eq. (|A59[) which scales nontriv- 
ially is the prefactor 6M/R cx R 2 as in the case of the 
pure top-hat. Eq. (|A43|) strictly holds for the compen- 
sated top-hat 



Ui = —M*i(8M,R) 



(A63) 



for the Newtonian and ip contributions. In particular, in 
the Vainshtein regime, the modification in U v from the 
pure top-hat is given by the constant piece in the profile 
which scales with e as in Fig. [T|>] 

In fact, this derivation unlike that for the pure top- 
hat is fully self-consistent in that the form of the profile, 
including the background contribution, is self-similar as 



the shells are removed. Furthermore the end result is a 
compensated top-hat of R — > 0, i.e. an unperturbed uni- 
verse with no source to ip. In this view the binding energy 
associated with SM is the energy required to eliminate 
the mass perturbation rather than the mass. 

Unfortunately, this definition still docs not fully resolve 
the ambiguities associated with the potential energy def- 
inition. We still need to account for the potential energy 




FIG. 16: Reduction in the surface potential due to compensa- 
tion in the Vainshtein limit, C E /Co, as a function of e = V*/8. 
In the limit e <C 1, y* <C y e and the compensation has lit- 
tle effect on the ip profile. This quantity also controls the 
reduction in binding energy. 



due to the background density. In particular, in this def- 
inition the Newtonian binding energy also only accounts 
for the perturbation and hence using it in the definition 
of total energy would not obey strict energy conservation 
during the collapse, especially for S < 1. We have seen 
that in the literature the definition of potential energy 
is mainly used in conjunction with energy conservation 
to simplify the calculation of the virial radius. In this 
context, the original definition of binding energy as that 
of a pure top-hat with no contribution from the exterior 
is more useful. In this case it is important to keep in 
mind however that modified forces as well as evolving 
dark energy density generally imply a violation of energy 
conservation (§ IA 5|) . 
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